home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / source / snip9503 / sunriset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-14  |  21.2 KB  |  509 lines
  1. /*
  2.  
  3. SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
  4.              the length of the day at any date and latitude
  5.  
  6. Written as DAYLEN.C, 1989-08-16
  7.  
  8. Modified to SUNRISET.C, 1992-12-01
  9.  
  10. (c) Paul Schlyter, 1989, 1992
  11.  
  12. Released to the public domain by Paul Schlyter, December 1992
  13.  
  14. */
  15.  
  16.  
  17. #include <stdio.h>
  18. #include <math.h>
  19.  
  20.  
  21. /* A macro to compute the number of days elapsed since 2000 Jan 0.0 */
  22. /* (which is equal to 1999 Dec 31, 0h UT)                           */
  23.  
  24. #define days_since_2000_Jan_0(y,m,d) \
  25.     (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
  26.  
  27. /* Some conversion factors between radians and degrees */
  28.  
  29. #ifndef PI
  30.  #define PI        3.1415926535897932384
  31. #endif
  32.  
  33. #define RADEG     ( 180.0 / PI )
  34. #define DEGRAD    ( PI / 180.0 )
  35.  
  36. /* The trigonometric functions in degrees */
  37.  
  38. #define sind(x)  sin((x)*DEGRAD)
  39. #define cosd(x)  cos((x)*DEGRAD)
  40. #define tand(x)  tan((x)*DEGRAD)
  41.  
  42. #define atand(x)    (RADEG*atan(x))
  43. #define asind(x)    (RADEG*asin(x))
  44. #define acosd(x)    (RADEG*acos(x))
  45. #define atan2d(y,x) (RADEG*atan2(y,x))
  46.  
  47.  
  48. /* Following are some macros around the "workhorse" function __daylen__ */
  49. /* They mainly fill in the desired values for the reference altitude    */
  50. /* below the horizon, and also selects whether this altitude should     */
  51. /* refer to the Sun's center or its upper limb.                         */
  52.  
  53.  
  54. /* This macro computes the length of the day, from sunrise to sunset. */
  55. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  56. /* 35 arc minutes below the horizon (this accounts for the refraction */
  57. /* of the Earth's atmosphere).                                        */
  58. #define day_length(year,month,day,lon,lat)  \
  59.         __daylen__( year, month, day, lon, lat, -35.0/60.0, 1 )
  60.  
  61. /* This macro computes the length of the day, including civil twilight. */
  62. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  63. /* the horizon.                                                         */
  64. #define day_civil_twilight_length(year,month,day,lon,lat)  \
  65.         __daylen__( year, month, day, lon, lat, -6.0, 0 )
  66.  
  67. /* This macro computes the length of the day, incl. nautical twilight.  */
  68. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  69. /* below the horizon.                                                   */
  70. #define day_nautical_twilight_length(year,month,day,lon,lat)  \
  71.         __daylen__( year, month, day, lon, lat, -12.0, 0 )
  72.  
  73. /* This macro computes the length of the day, incl. astronomical twilight. */
  74. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  75. /* below the horizon.                                                      */
  76. #define day_astronomical_twilight_length(year,month,day,lon,lat)  \
  77.         __daylen__( year, month, day, lon, lat, -18.0, 0 )
  78.  
  79.  
  80. /* This macro computes times for sunrise/sunset.                      */
  81. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  82. /* 35 arc minutes below the horizon (this accounts for the refraction */
  83. /* of the Earth's atmosphere).                                        */
  84. #define sun_rise_set(year,month,day,lon,lat,rise,set)  \
  85.         __sunriset__( year, month, day, lon, lat, -35.0/60.0, 1, rise, set )
  86.  
  87. /* This macro computes the start and end times of civil twilight.       */
  88. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  89. /* the horizon.                                                         */
  90. #define civil_twilight(year,month,day,lon,lat,start,end)  \
  91.         __sunriset__( year, month, day, lon, lat, -6.0, 0, start, end )
  92.  
  93. /* This macro computes the start and end times of nautical twilight.    */
  94. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  95. /* below the horizon.                                                   */
  96. #define nautical_twilight(year,month,day,lon,lat,start,end)  \
  97.         __sunriset__( year, month, day, lon, lat, -12.0, 0, start, end )
  98.  
  99. /* This macro computes the start and end times of astronomical twilight.   */
  100. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  101. /* below the horizon.                                                      */
  102. #define astronomical_twilight(year,month,day,lon,lat,start,end)  \
  103.         __sunriset__( year, month, day, lon, lat, -18.0, 0, start, end )
  104.  
  105.  
  106. /* Function prototypes */
  107.  
  108. double __daylen__( int year, int month, int day, double lon, double lat,
  109.                    double altit, int upper_limb );
  110.  
  111. int __sunriset__( int year, int month, int day, double lon, double lat,
  112.                   double altit, int upper_limb, double *rise, double *set );
  113.  
  114. void sunpos( double d, double *lon, double *r );
  115.  
  116. void sun_RA_dec( double d, double *RA, double *dec, double *r );
  117.  
  118. double revolution( double x );
  119.  
  120. double rev180( double x );
  121.  
  122. double GMST0( double d );
  123.  
  124.  
  125.  
  126. /* A small test program */
  127.  
  128. main()
  129. {
  130.       int year,month,day;
  131.       double lon, lat;
  132.       double daylen, civlen, nautlen, astrlen;
  133.       double rise, set, civ_start, civ_end, naut_start, naut_end,
  134.              astr_start, astr_end;
  135.       int    rs, civ, naut, astr;
  136.  
  137.       printf( "Longitude (+ is east) and latitude (+ is north) : " );
  138.       scanf( "%lf %lf", &lon, &lat );
  139.  
  140.       for(;;)
  141.       {
  142.             printf( "Input date ( yyyy mm dd ) (ctrl-C exits): " );
  143.             scanf( "%d %d %d", &year, &month, &day );
  144.  
  145.             daylen  = day_length(year,month,day,lon,lat);
  146.             civlen  = day_civil_twilight_length(year,month,day,lon,lat);
  147.             nautlen = day_nautical_twilight_length(year,month,day,lon,lat);
  148.             astrlen = day_astronomical_twilight_length(year,month,day,
  149.                   lon,lat);
  150.  
  151.             printf( "Day length:                 %5.2f hours\n", daylen );
  152.             printf( "With civil twilight         %5.2f hours\n", civlen );
  153.             printf( "With nautical twilight      %5.2f hours\n", nautlen );
  154.             printf( "With astronomical twilight  %5.2f hours\n", astrlen );
  155.             printf( "Length of twilight: civil   %5.2f hours\n",
  156.                   (civlen-daylen)/2.0);
  157.             printf( "                  nautical  %5.2f hours\n",
  158.                   (nautlen-daylen)/2.0);
  159.             printf( "              astronomical  %5.2f hours\n",
  160.                   (astrlen-daylen)/2.0);
  161.  
  162.             rs   = sun_rise_set         ( year, month, day, lon, lat,
  163.                                           &rise, &set );
  164.             civ  = civil_twilight       ( year, month, day, lon, lat,
  165.                                           &civ_start, &civ_end );
  166.             naut = nautical_twilight    ( year, month, day, lon, lat,
  167.                                           &naut_start, &naut_end );
  168.             astr = astronomical_twilight( year, month, day, lon, lat,
  169.                                           &astr_start, &astr_end );
  170.  
  171.             printf( "Sun at south %5.2fh UT\n", (rise+set)/2.0 );
  172.  
  173.             switch( rs )
  174.             {
  175.                 case 0:
  176.                     printf( "Sun rises %5.2fh UT, sets %5.2fh UT\n",
  177.                              rise, set );
  178.                     break;
  179.                 case +1:
  180.                     printf( "Sun above horizon\n" );
  181.                     break;
  182.                 case -1:
  183.                     printf( "Sun below horizon\n" );
  184.                     break;
  185.             }
  186.  
  187.             switch( civ )
  188.             {
  189.                 case 0:
  190.                     printf( "Civil twilight starts %5.2fh, "
  191.                             "ends %5.2fh UT\n", civ_start, civ_end );
  192.                     break;
  193.                 case +1:
  194.                     printf( "Never darker than civil twilight\n" );
  195.                     break;
  196.                 case -1:
  197.                     printf( "Never as bright as civil twilight\n" );
  198.                     break;
  199.             }
  200.  
  201.             switch( naut )
  202.             {
  203.                 case 0:
  204.                     printf( "Nautical twilight starts %5.2fh, "
  205.                             "ends %5.2fh UT\n", naut_start, naut_end );
  206.                     break;
  207.                 case +1:
  208.                     printf( "Never darker than nautical twilight\n" );
  209.                     break;
  210.                 case -1:
  211.                     printf( "Never as bright as nautical twilight\n" );
  212.                     break;
  213.             }
  214.  
  215.             switch( astr )
  216.             {
  217.                 case 0:
  218.                     printf( "Astronomical twilight starts %5.2fh, "
  219.                             "ends %5.2fh UT\n", astr_start, astr_end );
  220.                     break;
  221.                 case +1:
  222.                     printf( "Never darker than astronomical twilight\n" );
  223.                     break;
  224.                 case -1:
  225.                     printf( "Never as bright as astronomical twilight\n" );
  226.                     break;
  227.             }
  228.  
  229.       }
  230.       return 0;
  231. }
  232.  
  233.  
  234. /* The "workhorse" function for sun rise/set times */
  235.  
  236. int __sunriset__( int year, int month, int day, double lon, double lat,
  237.                   double altit, int upper_limb, double *trise, double *tset )
  238. /***************************************************************************/
  239. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  240. /*       Eastern longitude positive, Western longitude negative       */
  241. /*       Northern latitude positive, Southern latitude negative       */
  242. /*       The longitude value IS critical in this function!            */
  243. /*       altit = the altitude which the Sun should cross              */
  244. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  245. /*               for civil, -12 degrees for nautical and -18          */
  246. /*               degrees for astronomical twilight.                   */
  247. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  248. /*               Set to non-zero (e.g. 1) when computing rise/set     */
  249. /*               times, and to zero when computing start/end of       */
  250. /*               twilight.                                            */
  251. /*        *rise = where to store the rise time                        */
  252. /*        *set  = where to store the set  time                        */
  253. /*                Both times are relative to the specified altitude,  */
  254. /*                and thus this function can be used to comupte       */
  255. /*                various twilight times, as well as rise/set times   */
  256. /* Return value:  0 = sun rises/sets this day, times stored at        */
  257. /*                    *trise and *tset.                               */
  258. /*               +1 = sun above the specified "horizon" 24 hours.     */
  259. /*                    *trise set to time when the sun is at south,    */
  260. /*                    minus 12 hours while *tset is set to the south  */
  261. /*                    time plus 12 hours. "Day" length = 24 hours     */
  262. /*               -1 = sun is below the specified "horizon" 24 hours   */
  263. /*                    "Day" length = 0 hours, *trise and *tset are    */
  264. /*                    both set to the time when the sun is at south.  */
  265. /*                                                                    */
  266. /**********************************************************************/
  267. {
  268.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  269.       sr,         /* Solar distance, astronomical units */
  270.       sRA,        /* Sun's Right Ascension */
  271.       sdec,       /* Sun's declination */
  272.       sradius,    /* Sun's apparent radius */
  273.       t,          /* Diurnal arc */
  274.       tsouth,     /* Time when Sun is at south */
  275.       sidtime;    /* Local sidereal time */
  276.  
  277.       int rc = 0; /* Return cde from function - usually 0 */
  278.  
  279.       /* Compute d of 12h local mean solar time */
  280.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  281.  
  282.       /* Compute local sideral time of this moment */
  283.       sidtime = revolution( GMST0(d) + 180.0 + lon );
  284.  
  285.       /* Compute Sun's RA + Decl at this moment */
  286.       sun_RA_dec( d, &sRA, &sdec, &sr );
  287.  
  288.       /* Compute time when Sun is at south - in hours UT */
  289.       tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
  290.  
  291.       /* Compute the Sun's apparent radius, degrees */
  292.       sradius = 0.2666 / sr;
  293.  
  294.       /* Do correction to upper limb, if necessary */
  295.       if ( upper_limb )
  296.             altit -= sradius;
  297.  
  298.       /* Compute the diurnal arc that the Sun traverses to reach */
  299.       /* the specified altitide altit: */
  300.       {
  301.             double cost;
  302.             cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
  303.                   ( cosd(lat) * cosd(sdec) );
  304.             if ( cost >= 1.0 )
  305.                   rc = -1, t = 0.0;       /* Sun always below altit */
  306.             else if ( cost <= -1.0 )
  307.                   rc = +1, t = 12.0;      /* Sun always above altit */
  308.             else
  309.                   t = acosd(cost)/15.0;   /* The diurnal arc, hours */
  310.       }
  311.  
  312.       /* Store rise and set times - in hours UT */
  313.       *trise = tsouth - t;
  314.       *tset  = tsouth + t;
  315.  
  316.       return rc;
  317. }  /* __sunriset__ */
  318.  
  319.  
  320.  
  321. /* The "workhorse" function */
  322.  
  323.  
  324. double __daylen__( int year, int month, int day, double lon, double lat,
  325.                    double altit, int upper_limb )
  326. /**********************************************************************/
  327. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  328. /*       Eastern longitude positive, Western longitude negative       */
  329. /*       Northern latitude positive, Southern latitude negative       */
  330. /*       The longitude value is not critical. Set it to the correct   */
  331. /*       longitude if you're picky, otherwise set to to, say, 0.0     */
  332. /*       The latitude however IS critical - be sure to get it correct */
  333. /*       altit = the altitude which the Sun should cross              */
  334. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  335. /*               for civil, -12 degrees for nautical and -18          */
  336. /*               degrees for astronomical twilight.                   */
  337. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  338. /*               Set to non-zero (e.g. 1) when computing day length   */
  339. /*               and to zero when computing day+twilight length.      */
  340. /**********************************************************************/
  341. {
  342.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  343.       obl_ecl,    /* Obliquity (inclination) of Earth's axis */
  344.       sr,         /* Solar distance, astronomical units */
  345.       slon,       /* True solar longitude */
  346.       sin_sdecl,  /* Sine of Sun's declination */
  347.       cos_sdecl,  /* Cosine of Sun's declination */
  348.       sradius,    /* Sun's apparent radius */
  349.       t;          /* Diurnal arc */
  350.  
  351.       /* Compute d of 12h local mean solar time */
  352.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  353.  
  354.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  355.       obl_ecl = 23.4393 - 3.563E-7 * d;
  356.  
  357.       /* Compute Sun's position */
  358.       sunpos( d, &slon, &sr );
  359.  
  360.       /* Compute sine and cosine of Sun's declination */
  361.       sin_sdecl = sind(obl_ecl) * sind(slon);
  362.       cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
  363.  
  364.       /* Compute the Sun's apparent radius, degrees */
  365.       sradius = 0.2666 / sr;
  366.  
  367.       /* Do correction to upper limb, if necessary */
  368.       if ( upper_limb )
  369.             altit -= sradius;
  370.  
  371.       /* Compute the diurnal arc that the Sun traverses to reach */
  372.       /* the specified altitide altit: */
  373.       {
  374.             double cost;
  375.             cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
  376.                   ( cosd(lat) * cos_sdecl );
  377.             if ( cost >= 1.0 )
  378.                   t = 0.0;                      /* Sun always below altit */
  379.             else if ( cost <= -1.0 )
  380.                   t = 24.0;                     /* Sun always above altit */
  381.             else  t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
  382.       }
  383.       return t;
  384. }  /* __daylen__ */
  385.  
  386.  
  387. /* This function computes the Sun's position at any instant */
  388.  
  389. void sunpos( double d, double *lon, double *r )
  390. /******************************************************/
  391. /* Computes the Sun's ecliptic longitude and distance */
  392. /* at an instant given in d, number of days since     */
  393. /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
  394. /* computed, since it's always very near 0.           */
  395. /******************************************************/
  396. {
  397.       double M,         /* Mean anomaly of the Sun */
  398.              w,         /* Mean longitude of perihelion */
  399.                         /* Note: Sun's mean longitude = M + w */
  400.              e,         /* Eccentricity of Earth's orbit */
  401.              E,         /* Eccentric anomaly */
  402.              x, y,      /* x, y coordinates in orbit */
  403.              v;         /* True anomaly */
  404.  
  405.       /* Compute mean elements */
  406.       M = revolution( 356.0470 + 0.9856002585 * d );
  407.       w = 282.9404 + 4.70935E-5 * d;
  408.       e = 0.016709 - 1.151E-9 * d;
  409.  
  410.       /* Compute true longitude and radius vector */
  411.       E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
  412.             x = cosd(E) - e;
  413.       y = sqrt( 1.0 - e*e ) * sind(E);
  414.       *r = sqrt( x*x + y*y );              /* Solar distance */
  415.       v = atan2d( y, x );                  /* True anomaly */
  416.       *lon = v + w;                        /* True solar longitude */
  417.       if ( *lon >= 360.0 )
  418.             *lon -= 360.0;                   /* Make it 0..360 degrees */
  419. }
  420.  
  421. void sun_RA_dec( double d, double *RA, double *dec, double *r )
  422. {
  423.       double lon, obl_ecl, x, y, z;
  424.  
  425.       /* Compute Sun's ecliptical coordinates */
  426.       sunpos( d, &lon, r );
  427.  
  428.       /* Compute ecliptic rectangular coordinates (z=0) */
  429.       x = *r * cosd(lon);
  430.       y = *r * sind(lon);
  431.  
  432.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  433.       obl_ecl = 23.4393 - 3.563E-7 * d;
  434.  
  435.       /* Convert to equatorial rectangular coordinates - x is uchanged */
  436.       z = y * sind(obl_ecl);
  437.       y = y * cosd(obl_ecl);
  438.  
  439.       /* Convert to spherical coordinates */
  440.       *RA = atan2d( y, x );
  441.       *dec = atan2d( z, sqrt(x*x + y*y) );
  442.  
  443. }  /* sun_RA_dec */
  444.  
  445.  
  446. /******************************************************************/
  447. /* This function reduces any angle to within the first revolution */
  448. /* by subtracting or adding even multiples of 360.0 until the     */
  449. /* result is >= 0.0 and < 360.0                                   */
  450. /******************************************************************/
  451.  
  452. #define INV360    ( 1.0 / 360.0 )
  453.  
  454. double revolution( double x )
  455. /*****************************************/
  456. /* Reduce angle to within 0..360 degrees */
  457. /*****************************************/
  458. {
  459.       return( x - 360.0 * floor( x * INV360 ) );
  460. }  /* revolution */
  461.  
  462. double rev180( double x )
  463. /*********************************************/
  464. /* Reduce angle to within +180..+180 degrees */
  465. /*********************************************/
  466. {
  467.       return( x - 360.0 * floor( x * INV360 + 0.5 ) );
  468. }  /* revolution */
  469.  
  470.  
  471. /*******************************************************************/
  472. /* This function computes GMST0, the Greenwhich Mean Sidereal Time */
  473. /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
  474. /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
  475. /* time of the day.  I've generelized GMST0 as well, and define it */
  476. /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
  477. /* other times than 0h UT as well.  While this sounds somewhat     */
  478. /* contradictory, it is very practical:  instead of computing      */
  479. /* GMST like:                                                      */
  480. /*                                                                 */
  481. /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
  482. /*                                                                 */
  483. /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
  484. /* computes:                                                       */
  485. /*                                                                 */
  486. /*  GMST = GMST0 + UT                                              */
  487. /*                                                                 */
  488. /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
  489. /* Defined in this way, GMST0 will increase with about 4 min a     */
  490. /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
  491. /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
  492. /* (if we neglect aberration, which amounts to 20 seconds of arc   */
  493. /* or 1.33 seconds of time)                                        */
  494. /*                                                                 */
  495. /*******************************************************************/
  496.  
  497. double GMST0( double d )
  498. {
  499.       double sidtim0;
  500.       /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
  501.       /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
  502.       /* add these numbers, I'll let the C compiler do it for me.  */
  503.       /* Any decent C compiler will add the constants at compile   */
  504.       /* time, imposing no runtime or code overhead.               */
  505.       sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
  506.                           ( 0.9856002585 + 4.70935E-5 ) * d );
  507.       return sidtim0;
  508. }  /* GMST0 */
  509.